home *** CD-ROM | disk | FTP | other *** search
/ Eagles Nest BBS 8 / Eagles_Nest_Mac_Collection_Disc_8.TOAST / Developer Tools⁄Additions / MacScheme20 / Mathlib / re.scm < prev    next >
Encoding:
Text File  |  1989-04-27  |  8.1 KB  |  225 lines  |  [TEXT/????]
  1. ;;; $Header: re.scm,v 1.4 88/01/26 09:04:46 GMT gjs Exp $
  2. ;;; Richardson Extrapolation of a sequence, represented as a stream.
  3.  
  4. (if-mit
  5.  (declare (usual-integrations = + - * /
  6.                  zero? 1+ -1+
  7.                  ;; truncate round floor ceiling
  8.                  sqrt exp log sin cos)))
  9.  
  10. ;;; Assumption: The sequence is the values of a function at successive
  11. ;;; reciprocal powers of 2.  The function is assumed to be analytic in
  12. ;;; a region around h=0 with a power-series expansion.  The first
  13. ;;; non-constant term of the series has power p and successive terms
  14. ;;; are separated by powers of q.
  15.  
  16. (define (make-zeno-sequence R h)
  17.   (cons-stream (R h) (make-zeno-sequence R (/ h 2))))
  18.  
  19. (define (accelerate-zeno-sequence seq p)
  20.   (let* ((2^p (expt 2 p)) (2^p-1 (- 2^p 1)))
  21.     (map-streams (lambda (Rh Rh/2)
  22.                     (/ (- (* 2^p Rh/2) Rh)
  23.                        2^p-1))
  24.                  seq
  25.                  (tail seq))))
  26.  
  27. (define (make-zeno-tableau seq p q)
  28.   (define (sequences seq order)
  29.     (cons-stream seq
  30.              (sequences (accelerate-zeno-sequence seq order)
  31.                  (+ order q))))
  32.   (sequences seq p))
  33.  
  34. (define (first-terms-of-zeno-tableau tableau)
  35.   (map-stream head tableau))
  36.  
  37. (define (richardson-sequence seq p q)
  38.   (first-terms-of-zeno-tableau (make-zeno-tableau seq p q)))
  39.  
  40. (define (stream-limit s tolerance . optionals)
  41.   (let ((M (default-lookup 'MAXTERMS -1 optionals))
  42.         (ok? (default-lookup 'CONVERGENCE-TEST close-enuf? optionals)))
  43.     (let loop ((s s) (count 2))
  44.       (let* ((h1 (head s)) (t (tail s)) (h2 (head t)))
  45.         (cond ((ok? h1 h2 tolerance) h2)
  46.               ((and (positive? M) (>= count M))
  47.            ((default-lookup 'IFFAIL 
  48.                                 (lambda args 
  49.                                   (error "STREAM-LIMIT not reached" M))
  50.                                 optionals)
  51.             'stream-limit
  52.                 'not-converged
  53.         M
  54.         h1
  55.         h2))
  56.               (else (loop t (+ count 1))))))))
  57.  
  58. (define (richardson-limit f start-h ord inc tolerance . opts)
  59.   (apply stream-limit (richardson-sequence (make-zeno-sequence f start-h)
  60.                                            ord
  61.                                            inc)
  62.                       tolerance
  63.               opts))
  64.  
  65. (define ord-estimate-stream
  66.   (let ((log2 (log 2)))
  67.     (lambda (s)
  68.       (let* ((s0 (head s)) (st (tail s))
  69.              (s1 (head st)) (s2 (head (tail st))))
  70.     (cons-stream (/ (log (/ (- s0 s1) (- s1 s2))) log2)
  71.                  (ord-estimate-stream st))))))
  72.  
  73. (define (guess-integer-convergent s)
  74.   (let* ((s0 (head s)) (st (tail s)) (s1 (head st)) (s2 (head (tail st))))
  75.     (let* ((N (round s0))
  76.            (e2 (abs (- s2 N))) (e1 (abs (- s1 N))) (e0 (abs (- s0 N))))
  77.       (if (and (<= e2 e1) (<= e1 e0))
  78.           N
  79.           (guess-integer-convergent st)))))
  80.  
  81. (define (guess-ord-inc-etc s) ;tries to return { ord inc inc inc ... }
  82.   (let ((psequence (ord-estimate-stream s)))
  83.     (cons-stream (guess-integer-convergent psequence)
  84.                  (guess-ord-inc-etc psequence))))
  85.  
  86. (define (guess-ord-and-inc s)
  87.   (let ((g (guess-ord-inc-etc s)))
  88.     ;; ... and if we make it this far ...
  89.     (list (head g) (head (tail g)))))
  90.  
  91. ;;; For example, we can make a Richardson extrapolation to get the
  92. ;;; first derivative of a function to a required tolerance.
  93.         
  94. (define rderiv
  95.   (let ((log2 (log 2)))
  96.     (lambda (f tolerance)
  97.       (lambda (x)
  98.     (let* ((plausible-h (* 0.5 (abs x)))
  99.            (h (if (zero? plausible-h) 0.1 plausible-h))
  100.            (delta (- (f (+ x h)) (f (- x h))))
  101.            (roundoff (* *machine-epsilon*
  102.                 (+ 1 (floor (abs (/ (f x)
  103.                         (if (zero? delta)
  104.                             1
  105.                             delta)))))))
  106.            (n (floor (/ (log (/ tolerance roundoff)) log2))))
  107.       (richardson-limit (lambda (dx)
  108.                   (/ (- (f (+ x dx))
  109.                     (f (- x dx)))
  110.                  (* 2 dx)))
  111.                 h
  112.                 2
  113.                 2
  114.                 tolerance
  115.                 'MAXTERMS (+ n 1)))))))
  116.  
  117. (define rsecond-deriv
  118.   (let ((log4 (log 4)))
  119.     (lambda (f tolerance)
  120.       (lambda (x)
  121.     (let* ((plausible-h (* 0.5 (abs x)))
  122.            (h (if (zero? plausible-h) 0.1 plausible-h))
  123.                (2fx (* 2 (f x)))
  124.            (delta (+ (f (+ x h h)) (- 2fx) (f (- x h h))))
  125.            (roundoff (* *machine-epsilon*
  126.                 (+ 1 (floor (abs (/ 2fx
  127.                         (if (zero? delta)
  128.                             1
  129.                             delta)))))))
  130.            (n (floor (/ (log (/ tolerance roundoff)) log4))))
  131.       (richardson-limit (lambda (h)
  132.                   (/ (+ (f (+ x h h))
  133.                     (* -2 (f x))
  134.                     (f (- x h h)))
  135.                  (* 4 h h)))
  136.                 h
  137.                 2
  138.                 2
  139.                 tolerance
  140.                 'MAXTERMS (+ n 1)))))))
  141.  
  142. ;;; Romberg integration is the result of a Richardson extrapolation
  143. ;;; of the trapezoid method.
  144.  
  145. (define (rinteg f a b tolerance)
  146.   (define (trapezoid-sums f a b)
  147.     (define (next-S S n)
  148.       (let* ((h (/ (- b a) 2 n))
  149.          (fx (lambda(i) (f (+ a (* (+ i i -1) h))))))
  150.     (+ (/ S 2) (* h (sigma fx 1 n)))))
  151.     (define (S-and-n-stream S n)
  152.       (cons-stream (list S n)
  153.            (S-and-n-stream (next-S S n) (* n 2))))
  154.     (let* ((h (- b a))
  155.        (S (* (/ h 2) (+ (f a) (f b)))))
  156.       (map-stream car (S-and-n-stream S 1))))
  157.   (stream-limit (richardson-sequence (trapezoid-sums f a b) 
  158.                                      2
  159.                          2)
  160.                 tolerance))
  161.  
  162. ;;; Given a sequence (stream) of abscissas and ordinates for a function,
  163. ;;; the following procedure constructs a sequence of constant terms of 
  164. ;;; polynomials that interpolate the successive initial segments of the 
  165. ;;; given sequences.
  166. ;;; This is done by incrementally constructing a Neville-like tableau
  167. ;;; for the interpolating polynomials, specialized for argument 0.
  168.  
  169. (define polynomial-extrapolation
  170.   (letrec ((pt-lp
  171.              (lambda (xs ys xstate ystate)
  172.                (let ((x (head xs)))
  173.                  (letrec ((poly-lp
  174.                             (lambda (y xstate ystate)
  175.                               (if (null? ystate)
  176.                                   (list y)
  177.                                   (cons y
  178.                                         (poly-lp (/ (- (* y (car xstate))
  179.                                                        (* x (car ystate)))
  180.                                                     (- (car xstate) x))
  181.                                                  (cdr xstate)
  182.                                                  (cdr ystate)))))))
  183.                    (let ((new (poly-lp (head ys) xstate ystate)))
  184.                      (cons-stream (car (last-pair new))
  185.                                   (pt-lp (tail xs)
  186.                                          (tail ys)
  187.                                          (cons x xstate)
  188.                                          new))))))))
  189.     (lambda (abscissas ordinates)
  190.       (pt-lp abscissas ordinates '() '()))))
  191.  
  192. ;;; Given a sequence (stream) of abscissas and ordinates for a function,
  193. ;;; the following procedure constructs a sequence of constant terms of 
  194. ;;; rational functions that interpolate the successive initial segments
  195. ;;; of the given sequences.
  196. ;;; This is done by incrementally constructing a Bulirsch-Stoer tableau
  197. ;;; for the interpolating rational functions, specialized for argument 0.
  198.  
  199. (define rational-extrapolation
  200.   (letrec ((pt-lp
  201.              (lambda (xs ys xstate ystate)
  202.                (let ((x (head xs)))
  203.                  (letrec ((rat-lp
  204.                             (lambda (y xstate ystate z)
  205.                               (if (null? ystate)
  206.                                   (list y)
  207.                   (let ((u (* (/ (car xstate) x) 
  208.                                               (- (car ystate) z))))
  209.                                     (cons y
  210.                                       (rat-lp
  211.                                         (/ (+ (* y (- u (car ystate)))
  212.                                               (* z (car ystate)))
  213.                                            (- u (- y z)))
  214.                                         (cdr xstate)
  215.                                         (cdr ystate)
  216.                     (car ystate))))))))
  217.                    (let ((new (rat-lp (head ys) xstate ystate 0)))
  218.                      (cons-stream (car (last-pair new))
  219.                                   (pt-lp (tail xs)
  220.                                          (tail ys)
  221.                                          (cons x xstate)
  222.                                          new))))))))
  223.     (lambda (abscissas ordinates)
  224.       (pt-lp abscissas ordinates '() '()))))
  225.